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AN  IMPLICIT  SEMIANALYTIC  NUMERICAL  METHOD  FOR  THE 
SOLUTION  OF  NONEQUILIBRIUM  CHEMISTRY  PROBLEMS 

By  R.  A.  Graves,  Jr.,  P.  A.  Gnoffo,  and  R.  E.  Boughner 


INTRODUCTION 


Many  physical  phenomena  are  modeled  by  systems  of  linear  and/or  nonlinear 
ordinary  differential  equations  (see  for  exanqple  references  1  to  5)  which  are 
defined  as  stiff  systems  when  a  large  spread  in  negative  eigenvalues  exists. 

Such  stiff  systems  commonly  arise  in  nonequilibrium  chemistry  problems  involving 
kinetic  and  photochemical  reactions,  "nie  governing  equations  for  these  stiff 
systems  are  difficult  to  solve  numerically  using  classical  techniques  because 
the  error  growth  is  rapid,  and  unless  the  equations  are  integrated  using  a  very 
small  time  step,  the  results  can  be  meaningless.  To  alleviate  the  problems 
involved  with  stiff  systems,  a  great  deal  of  effort  has  been  expended  in 
developing  numerical  solution  techniques,  both  explicit  and  implicit,  for 
stiff  ordinary  differential  equations.  References  6  to  9  review  some  of  the 
more  popular  numerical  methods  and  present  the  results  of  numerical  comparisons 
between  the  methods.  A  generalized  conclusion  resulting  from  the  studies  of 
references  6  to  8  is  that  the  implicit  methods  are  more  desirable  because  of 
their  increased  stability  and,  in  some  instances,  significantly  fewer  mathe¬ 
matical  operations.  In  these  and  other  studies,  a  rather  simple  (yet 
fundamental)  implicit  technique  was  not  investigated  because  these 
studies  used  a  generalized  equation  which  did  not  take  advantage  of 
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the  form  of  the  governing  conservation  equation  for  chemical  species.  The 
governing  conservation  equations  for  systems  of  chemical  reactions  can  generally 
be  written  in  the  form  of  first-order  ordinary  differential  equations.  These 
equations  can  be  solved  by  a  simple  implicit  semianalytic  technique  which  is 
derived  from  a  quadrature  solution  of  the  governing  equations.  This  method 
is  mathematically  simpler  than  most  implicit  methods  and  has  the  exponen¬ 
tial  nature  of  the  problem  embedded  in  the  solution. 

The  objective  of  this  paper  is  to  present  the  development  of  the  semi¬ 
analytic  technique  (SAT)  and  to  compare  its  efficiency  to  that  of  several  of 
the  more  popular  methods  available. 

SYMBOLS 

a,b  general  coefficients,  see  equation  (1) 

Cj,C2  curve  fit  coefficients,  see  equation  (5) 

^1*^2  curve  fit  coefficients,  see  equation  (6) 

HS  Hermite-Simpsoii 

Rn  [Yn  -  Y(Xjj)]/Y(Xn) 

h  step  size 

RK4  fourth-order  Runge-Kutta 

DEQ  Adams'  fourth-order  P-C 

TM  Treanor's  method 

DIFSYS  modified  midpoint  rule 

TR  Trapezoidal  rule 

TR-EX  Trapezoidal  rule  with  extrapolation 

CAL  Calahan's  method 

LNl  Liniger-Willoughby  -  Method  1 
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Liniger-Willoughby  -  Method  3 
semianalytic  technique 
eigenvalue 
calculated  value 
'•exact  value" 

transformed  coordinate,  see  equation  3a 

MATHEMATICAL  DEVELOPMENT 

The  governing  equation  for  the  conservation  of  chemical  species  in 
nonequilibrium  chemically  reacting  systems  can  generally  be  written  in  the 
form  of  a  first-order  ordinary  differential  equation  (see  ref.  10): 

dY. 

+  a(t)Yj^  =  b(t)  (1) 

where  a(t)  and  b(t)  generally  represent  the  loss  and  production  rates  of 
species  i,  respectively.  The  solution  of  equation  1  in  terms  of  quadratures 
is:  (This  procedure  is  similar  to  that  used  in  ref.  11.) 

tJ+1  t-^+l 

a(t)dt  ^J+1  -/  a(t')dt* 

+  /  b(t)e  ^  dt  (2) 

t'J 

This  equation  can  be  further  simplified  by  introducing  the  transformation 

tJ+1 

5  »  /  a(t')dt' 
dC  “  -a(t)dt 


J+1  J  tJ 
Y.  *  Y.  e  ^ 

1  1 
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hence,  equation  2  becomes: 
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(3) 


where 


/  a(t)dt 
td 


(4) 


For  the  least  complicated  case,  the  coefficients  a(t)  and  b(C)/a(C) 
can  be  approximated  by  linear  functions. 


a(t)  -  Cj  +  C2  (t  -  t*^) 


(5) 


where 

Cj  =  »(tJ) 

c  -  -  a(tJ) 

2  At 


bm  _ 

rtf)  ■ 


+  C2C 


(6) 


where 


S 


b(0) 

a(0) 


b(Ci)  _  b(0) 

C  =  a(0) 

2  AC 


It  should  be  noted  that  due  to  the  transformation  from  t  to  C  that 


b(Ci)  .  b(tJ) 
a(Ci)  a(td) 


Introducing  equation  6  into  equation  3  results  in: 


-  W  e“^^  +  (C,  +  C25)e"^dC 

X  X  1  X 


(7) 
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This  equation  can  now  be  integrated  by  parts  to  obtain  the  following  semi- 
analytic  implicit  result: 

T+i  T 

yp  =  YT  e  %  (C^  +  C2)(l  -  e  ")  -  (8) 

(Note:  To  have  stability  and  accuracy,  it  is  necessary  that  <  1.) 

Equation  8  is  semianalytic  in  nature  and  includes  the  inherent  exponential 
behavior  of  the  stiff  problem  directly  in  the  solution.  Equation  8  must  be 
solved  implicitly  (iteratively)  as  the  constants  and  depend  on  the 

conditions  at  the  advanced  time 

An  error  analysis  was  performed,  using  the  method  of  chapter  2  of  refer¬ 
ence  12,  to  deteirmine  the  errors  incurred  in  making  the  linear  approximations 
for  the  coefficients  a(t)  and  b(5)/a(C).  The  lowest  order  error  terms  are: 

E  =  -  ^  {A5V'(a>)  +  At^3»(n)(^l  ^  ^i^} 

where  y"iui)  is  the  second  derivative  of  the  ratio  b(C)/a(5)  on  the  interval 
0  <  0)  <  arid  6''(n)  is  the  second  derivative  of  a(t)  on  the  interval 

t*^  <  n  < 

Numerical  Experiments: 

System  I  (ref.  8) 

=  -0.04Yj  +  J0'^Y2Y3 

Y^  =  0.04Yj  -  10^Y2Y3  -  3  x  10^  Y^ 

Y^  =  3  X  10^  Y2 

Yj(0)  =  1  Y2(0)  =  0  Y3(0)  =  0 

This  system  is  nonlinear, and  no  exact  solution  for  this  system  exists. 

The  eigenvalues,  determined  from  the  Jacobian  matrix  of  the  system  at  X  =  0 
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are  X,  »  0,  X-  =  0  and  X_  =  -0.04.  X  changes  from  0.04  to  2405  for 

0  5  X  <  0.02.  The  eigenvalues  for  0  <  X  <  40  are  all  strictly  negative  or 
zero,  with  X^  *  0,  X2  ~  -10  ^  and  X^  -  -10^  to  -10^.  The  sharp  increase  in 
the  magnitude  of  X,  makes  this  a  particularly  difficult  stiff  system  to 
work  with.  In  addition  this  system  presents  some  starting  problems  for  SAT 
since  a2(0)  =  0  and  hence  meaningless.  To  circumvent  this 

problem  for  two  techniques  using  constant  h  were  tried:  first,  the 

Hermite-Simpson  method,  reference  13,  and  secondly,  the  Runge-Kutta  fourth-order 
method.  The  RK4  start  gave  the  best  results.  Table  I  gives  the  results  for 
this  system  on  the  CDC  6600  as  well  as  the  results  of  Lapidus  and  Seinfeld, 
reference  8,  for  the  IBM  7094.  (The  CDC  6600  is  approximately  10  times  faster 
than  the  IBM  7094.)  It  should  be  noted  that  due  to  the  nature  of  this  system 
YjCX)  was  calculated  by  Y^CX)  =  1  -  Y2(X)  -  Y^(X). 

The  most  successful  application  of  the  SAT  was  to  use  the  RK4  one  step 
to  obtain  Y^ (0.0005),  Y2 (0.0005),  and  Yj (0.0005)  and  then  use  the  SAT  with 

a  step  size  of  0.2  to  compute  the  solution  from  5  X  lO"^  <  X  <  40.  As  can  be 
seen  in  figures  1  and  2  the  linear  approximation  for  b(5)/a(C)  is  very 
accurate  for  this  system. 

System  II  (ref.  8) 

Y'  =  -200  (Y  -  F(X))  +  F'(X) 

Y(0)  =  10 

F(X)  =  10  -  (10  +  X)e'^ 

Exact  Solution  Y(X)  =  F(X)  +  lOe"^®®^ 

Here,  F(X)  is  a  slowly  decaying  solution  component  and  lOe"  decays 

rapidly.  The  large  negative  eigenvalue  of  -200  makes  exp(-200X)  negligible 
coB^ared  to  exp  (-X)  in  the  F(X)  component.  Results  using  SAT  on  the 
CDC  6600  are  compared  to  results  obtained  by  Lapidus  and  Seinfeld  in  Table  2. 
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At  equivalent  step  sizes,  SAT  produced  less  than  or  equal  to  the  error 

encountered  using  other  methods,  and  worked  faster  than  any  other  method  (times 
for  this  system  are  based  on  computing  over  a  range  0  <  X  <  15) 

System  III  (ref.  8) 

YJ  =  -O.lYj  -  49.9Y2 


Y*  = 

2 

-5OY2 

Y*  = 

3 

7OY2  -  I2OY3 

Yi(0) 

=  2  Y2(0)  =  : 

Exact 

Solution: 

Y^CX) 

-O.IX  -SOX 

=  e  +  e 

Y2(X) 

-SOX 
=  e 

YjCX) 

-SOX  -120X 

a  e  +  e 

Eigenvalues  Xj  =  -120 

YjCO)  =  2 


,  =  -50,  Xj  =  -0.1 

Because  the  solution  components  due  to  and  X2  decay  rapidly,  a 

stiff  method  which  was  not  restricted  by  the  magnitude  of  these  values  is 
desired.  The  SAT  (which  yields  the  exact  solution  of  Y2(X)  for  any  h  since 
the  linear  approximation  to  a(X)  and  b/a(5)  gives  the  true  variation 
of  these  functions)  is  compared  to  results  obtained  by  Lapidus  and  Seinfeld  in 
Table  3.  An  h  *  0.01  produced  results  which  were  better  than  the  results 
obtained  by  any  other  method.  However  as  the  step  size  increased,  the 
accuracy  dropped  off  rapidly  and  at  an  h  =  0.2,  the  solution  was  very  different 
from  the  exact  solution  (except,  of  course,  for  Y2(X)  which  remained  exact). 
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After  examining  the  problem,  it  was  found  that  on  an  interval 
of  0  -  X  -  0.2,  with  h  =  0.2,  the  linear  approximation  to  b/a(C)  was 
very  poor.  For  example,  b/a(X)  =  -400e  and  b^/a^CX)  =  (70/120)e 

(see  fig.  3).  The  effect  of  these  terras  on  b/a(C)  decays  rapidly  after 
X  *  0.2,  and  they  can  be  approximated  by  a  linear  function,  so  a  method  was 
tried  using  h  »  0.01  to  arrive  at  X  =  0.2  and  then  proceed  from  X  =  0.2 
with  h  =  0.2.  This  method  was  the  fastest  and  yielded  reasonable  results. 

System  IV 

yJ  =  O.8Y2  -  O.OlYj  -  lO^jY^Yj  +  lOYjYj  -  100YjY2 

y1  =  -0.8Y.,  -  10Y,Y-  +  loS.Y.  +1o'’y-Y. 

2  2  13  24  14 

Yl  =  O.OIY?  +  1oVy-Y-  +  2000Y^  -  lo'^Y.Y. 

3  1  123  4  14 

Y^  »  -I0S2Y4  "  20000Y^ 

Y^(0)  =  0.9;  Y2(0)  =  0.05;  Yj(0)  =  0.05;  Y^(0)  =  0 

No  exact  solution  for  this  nonlinear  system  was  obtained.  The  eigenvalues 

for  this  system,  calculated  from  the  Jacobian  using  values  of  Y  from  RK4, 

are  widely  separated  in  magnitudes.  All  of  the  eigenvalues  are  negative  or 

*0  on  a  range  6.146  x  10~^  <  X  <  7.36  x  10”^.  Typical  values  on  the  range  are 

at  X  =  7  X  10"^,  Xj  =  -1.017  x  10^,  X2  =  -4.979  x  lo"^,  X^  =  -2.7102  x  10^  and 

_g 

X^  =  -2.515  X  10  .  Results  for  this  system  appear  in  Table  4.  This  system 

was  only  stiff  for  a  short  time  and  none  of  the  methods  had  problems  with 
stability  on  a  range  0  <  X  <  2  x  10~5.  Using  RK4  with  h  =  1  X  10~®  as  a 
standard  of  comparison,  the  tables  indicate  that  for  any  given  step  size,  RK4 
was  more  accurate  than  any  implicit  method.  SAT  gave  accuracy  comparable  to 
other  implicit  methods  and  ran  at  approximately  the  same  speeds  as  these 
methods  for  equivalent  step  sizes.  In  this  system,  SAT  showed  no  advantage 


over  any  other  method.  However,  this  system  does  indicate  that  SAT  gives 
reasonable  results  for  nonstiff  nonlinear  systems. 

CONCLUDING  REMARKS 

As  developed  herein  the  crucial  approximation  is  the  linearization  of 
the  coefficients  within  the  integral  in  the  quadrature  form  which  allows  the 
semianalytic  form  to  be  obtained.  In  some  cases  this  approximation  is  very 
good  but  in  some  applications,  the  linear  approximation  can  be  in  error  for 
what  appear  to  be  not  unreasonable  time  steps.  As  demonstrated,  this  problem 
was  overcome  by  having  a  variable  time  step  which  is  small  in  the  region  where 
the  Linear  Approximation  is  in  error. 

Additionally,  quadratic  and  exponential  curve  fits  were  tried.  However 
these  approximations  produced  results  which  were  approximately  equivalent  to 
those  obtained  with  the  linear  approximation.  Because  the  linear  approxima¬ 
tion  was  the  simplest  to  program  and  because  of  the  consistently  good  results 
it  yielded,  it  was  chosen  as  the  most  desirable  approximation  evaluated. 

An  important  feature  of  the  semianalytic  technique  is  that  it  will  allow 
the  computation  of  nonequilibrium  chemical  systems  to  and  including  the  equi¬ 
librium  state.  For  systems  where  the  rates  are  large  (typical  of  approaching 
equilibrium)  the  SAT  equilibrium  condition  is  the  exact  solution  for 
equilibrium. 

As  demonstrated  in  the  example  problems,  the  semianalytic  technique  is 
both  rapid  and  accurate  and  should  be  applicable  to  those  stiff  problems  which 
can  be  modeled  by  an  equation  like  that  used  in  this  development. 
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TABLE  I.-  COMPARISON  OF  RESULTS  FOR  SYSTEM  I 


II 


TABLE  2.-  COMPARISON  OF  RESULTS  FOR  SYSTEM  II 
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(a)  Unstable 


TABLE  3.-  COMPARISON  OF  RESULTS  FOR  SYSTEM  III 
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TABLE  4.-  COMPARISON  OF  RESULTS  FOR  SYSTEM  IV 
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b)  Reaction  2 


Figure  !■  Variation  of  (a/b)  rjr  System  I  in  the  interval  .2. 


Linear  AppRoximiicjN  /  True  Variation 


Variation  of  (b/a)  for  System  III. 


